function groups = user_grouping(angles)
% USER_GROUPING 用户分组算法（Algorithm 1）
%   输入：
%     angles: 用户角度结构体数组（generate_user_angles输出）
%   输出：
%     groups: [K×1]向量，表示每个用户所属的组

    % ========== 1. 初始化参数 ==========
    K = length(angles);          % 用户总数
    P = length(angles(1).gamma_x); % 每条路径数
    
    % ========== 2. 构建相似度矩阵J（公式5） ==========
    J = zeros(K, K); % 初始化相似度矩阵
    
    for i = 1:K
        for j = i+1:K % 利用对称性减少计算量 因为 ρ12 = ρ21 的
            % 获取用户i和j的角度向量
            gamma_x_i = angles(i).gamma_x;
            gamma_y_i = angles(i).gamma_y;
            gamma_x_j = angles(j).gamma_x;
            gamma_y_j = angles(j).gamma_y;
            
            % 计算相似度（公式5）
            sim = 0;
            for p = 1:P
                for q = 1:P
                    dx = gamma_x_i(p) - gamma_x_j(q);
                    dy = gamma_y_i(p) - gamma_y_j(q);
                    sim = sim + exp(-(dx^2 + dy^2));
                end
            end
            
            % 对称赋值
            J(i, j) = sim;
            J(j, i) = sim;
        end
    end
    
    % ========== 3. 构建拉普拉斯矩阵L（公式6） ==========
    % 计算度向量τ_i = ∑_{j=1}^K J(i,j)
    tau = zeros(1,K);
    for i = 1 : K
        tau(i) = sum(J(i,:)) - J(i,i);
    end
%     tau = sum(J, 2);
%     
    % 构建对角矩阵N = diag(τ)
    N = diag(tau);
    
    % 拉普拉斯矩阵L = N - J
    L = N - J;
    
    % ========== 4. 特征分解确定组数G ==========
    % 特征分解（文献公式7）
    [V, D] = eig(L);
    
    % 容差阈值（处理数值误差）
    zero_eig_tol = 5;
    
    % 找到特征值接近0的索引
    zero_eig_idx = find((diag(D)) < zero_eig_tol);
    
    % 确定组数G（特征值0的重数）
    G = length(zero_eig_idx);
    
    % 处理特殊情况：无零特征值
    if G == 0
        warning('No zero eigenvalue found. Using single group.');
        G = 1;
        zero_eig_idx = 1; % 取最小特征值
    end
    
    % 提取特征向量矩阵Q_G（文献公式8）
    Q_G = V(:, zero_eig_idx);
    
    % ========== 5. K-means聚类分组 ==========
    % 应用K-means对Q_G的行向量聚类
    [groups, ~] = kmeans(Q_G, G, 'Replicates', 10, 'MaxIter', 1000);
    
    % ========== 6. 验证分组结果 ==========
%     validate_grouping(angles, groups, G);
end

%% ========== 子函数：分组结果验证 ==========
function validate_grouping(angles, groups, G)
    K = length(groups);
    
    % 验证组数
    if length(unique(groups)) ~= G
        error('Number of groups does not match eigenvalue multiplicity');
    end
    
    % 验证组内用户角度相似性
    for g = 1:G
        group_users = find(groups == g);
        
        % 计算组内角度范围
        all_azimuth = [];
        all_elevation = [];
        for u = group_users
            all_azimuth = [all_azimuth; angles(u).azimuth_angles];
            all_elevation = [all_elevation; angles(u).elevation_angles];
        end
        
        % 验证角度范围连续性（文献图6）
        azimuth_range = range(all_azimuth);
        elevation_range = range(all_elevation);
        
        if (azimuth_range > 30 || elevation_range > 25)
            warning('Group %d has large angle spread (Az: %.1f°, El: %.1f°)',...
                g, azimuth_range, elevation_range);
        end
    end
end